﻿/* Copyright 2007-2009 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 * 
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
using System;
using dnAnalytics.Statistics.Distributions;


namespace dnAnalytics.Mcmc
{
    /// <summary>
    /// Metropolis-Hastings sampling produces samples from distribition P by sampling from a proposal distribution Q
    /// and accepting/rejecting based on the density of P. Metropolis-Hastings sampling doesn't require that the
    /// proposal distribution Q is symmetric in comparison to <seealso cref="MetropolisSampler{T}"/>. It does need to
    /// be able to evaluate the proposal sampler's log density though. All densities are required to be in log space.
    /// 
    /// The Metropolis-Hastings sampler is a stateful sampler. It keeps track of where it currently is in the domain
    /// of the distribution P.
    /// </summary>
    /// <typeparam name="T">The type of samples this sampler produces.</typeparam>
    public class MetropolisHastingsSampler<T> : McmcSampler<T>
    {
        /// <summary>
        /// Evaluates the log density function of the target distribution.
        /// </summary>
        private readonly DensityLn<T> mPdfLnP;

        /// <summary>
        /// Evaluates the log transition probability for the proposal distribution.
        /// </summary>
        private readonly TransitionKernelLn<T> mKrnlQ;

        /// <summary>
        /// A function which samples from a proposal distribution.
        /// </summary>
        private readonly LocalProposalSampler<T> mProposal;

        /// <summary>
        /// The current location of the sampler.
        /// </summary>
        private T mCurrent;

        /// <summary>
        /// The log density at the current location.
        /// </summary>
        private double mCurrentDensityLn;

        /// <summary>
        /// The number of burn iterations between two samples.
        /// </summary>
        private int mBurnInterval;

        /// <summary>
        /// Constructs a new Metropolis-Hastings sampler using the default <see cref="System.Random"/> random 
        /// number generator. The burn interval will be set to 0.
        /// </summary>
        /// <param name="x0">The initial sample.</param>
        /// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
        /// <param name="krnlQ">The log transition probability for the proposal distribution.</param>
        /// <param name="proposal">A method that samples from the proposal distribution.</param>
        public MetropolisHastingsSampler(T x0, DensityLn<T> pdfLnP, TransitionKernelLn<T> krnlQ, LocalProposalSampler<T> proposal) :
            this(x0, pdfLnP, krnlQ, proposal, 0)
        {
        }

        /// <summary>
        /// Constructs a new Metropolis-Hastings sampler using the default <see cref="System.Random"/> random number generator. This
        /// constructor will set the burn interval.
        /// </summary>
        /// <param name="x0">The initial sample.</param>
        /// <param name="pdfLnP">The log density of the distribution we want to sample from.</param>
        /// <param name="krnlQ">The log transition probability for the proposal distribution.</param>
        /// <param name="proposal">A method that samples from the proposal distribution.</param>
        /// <param name="burnInterval">The number of iterations in between returning samples.</param>
        /// <exception cref="ArgumentOutOfRangeException">When the number of burnInterval iteration is negative.</exception>
        public MetropolisHastingsSampler(T x0, DensityLn<T> pdfLnP, TransitionKernelLn<T> krnlQ, LocalProposalSampler<T> proposal, int burnInterval)
        {
            mCurrent = x0;
            mCurrentDensityLn = pdfLnP(x0);
            mPdfLnP = pdfLnP;
            mKrnlQ = krnlQ;
            mProposal = proposal;
            BurnInterval = burnInterval;

            Burn(BurnInterval);
        }

        /// <summary>
        /// Gets or sets the number of iterations in between returning samples.
        /// </summary>
        /// <exception cref="ArgumentOutOfRangeException">When burn interval is negative.</exception>
        public int BurnInterval
        {
            get { return mBurnInterval; }

            set
            {
                if (value < 0)
                {
                    throw new ArgumentOutOfRangeException(Properties.Resources.ParameterCannotBeNegative);
                }
                mBurnInterval = value;
            }
        }

        /// <summary>
        /// This method runs the sampler for a number of iterations without returning a sample
        /// </summary>
        private void Burn(int n)
        {
            for (int i = 0; i < n; i++)
            {
                // Get a sample from the proposal.
                T next = mProposal(mCurrent);
                // Evaluate the density at the next sample.
                double p = mPdfLnP(next);
                // Evaluate the forward transition probability.
                double fwd = mKrnlQ(next, mCurrent);
                // Evaluate the backward transition probability
                double bwd = mKrnlQ(mCurrent, next);

                mSamples++;

                double acc = System.Math.Min(0.0, p + bwd - mCurrentDensityLn - fwd);
                if (acc == 0.0)
                {
                    mCurrent = next;
                    mCurrentDensityLn = p;
                    mAccepts++;
                }
                else if (Bernoulli.Sample(RandomNumberGenerator, System.Math.Exp(acc)) == 1)
                {
                    mCurrent = next;
                    mCurrentDensityLn = p;
                    mAccepts++;
                }
            }
        }

        /// <summary>
        /// Returns a sample from the distribution P.
        /// </summary>
        public override T Sample()
        {
            Burn(BurnInterval + 1);

            return mCurrent;
        }
    }
}